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We introduce a theoretical formalism to describe disorder-induced extrinsic scattering in slow- 
light photonic crystal waveguides. This work details and extends the optical scattering theory used 
in a recent Physical Review Letter [M. Patterson et al., Phys. Rev. Lett. 102, 103901 (2009)] to 
describe coherent scattering phenomena and successfully explain complex experimental measure- 
ments. Our presented theory, that combines Green function and coupled mode methods, allows 
one to self-consistently account for arbitrary multiple scattering for the propagating electric field 
and recover experimental features such as resonances near the band edge. The technique is fully 
three-dimensional and can calculate the effects of disorder on the propagating field over thousands 
of unit cells. As an application of this theory, we explore various sample lengths and disordered 
instances, and demonstrate the profound effect of multiple scattering in the waveguide transmission. 
The spectra yield rich features associated with disorder-induced localization and multiple scattering, 
which are shown to be exasperated in the slow light propagation regime. 
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Photonic crystal (PC) waveguides are structures 
formed by a line defect in an otherwise nominally per- 
fect photonic crystal lattice. PC slab waveguides are 
of particular interest because they can be fabricated us- 
ing high quality etching and lithography techniques. By 
guiding light using the photonic band gap of the sur- 
rounding crystal, strong transverse confinement on the 
order of a wavelength can be achieved. PC waveguides 
often exhibit a region of slow light propagation 1 '- which 
has potential applications as an optical delay line ,! or for 
enhanced light-matter interactions. 

It is now widely accepted that slow light propaga- 
tion enhances scattering from structural imperfections 
or fabrication disorder, leading to significant propaga- 
tion losses 4 ' 5 . Incoherent scattering theories that cal- 
culate the loss in a single waveguide period averaged 
over many nominally identical samples have predicted 
backscattcring and radiative loss to scale with the group 
velocity v g , as v~ 2 and v" 1 respectively 1 ' -8 . These ap- 
proximate loss-scaling relations have been confirmed ex- 
perimentally, e.g. 4 ' 9,10 , but they break down at low group 
velocities where multiple disorder-induced scattering be- 
comes significant. The simple scaling trends expected 
also typically do not include effects such as variation of 
the Bloch mode with wave vector or extrapolating the 
unit-cell loss to multiple waveguide periods, though re- 
cent work has included such effects within an incoherent 
scattering approach and shown a dramatic impact on the 
loss versus group velocity scaling rules 11 . Enhanced scat- 
tering losses in other material systems also occur in the 
slow light regime, for example, massive losses also occur 
in slow- light metamaterial waveguides 1- . 

In a recent Physical Review Letter 1 by Patterson 
et al., we extended previous theoretical incoherent- 
scattering work 1 ' to model coherent scattering over the 
entire length of a disordered waveguide instance, as 
schematically illustrated in Figure 1. This theory ex- 
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FIG. 1: Schematic of the problem under consideration: a 
plan view of a PC waveguide with scattering sites (stars) 
whose strength is enhanced by the vanishing group veloc- 
ity. Light injected from the left, indicated by the arrows, 
undergoes scattering at each of the sites leading to a com- 
plex interplay of forward and backward propagating waves. 
Ultimately, most of the light is back scattered and the trans- 
mission is low. In the full calculation, the scatterers are con- 
tinuously distributed throughout the system and we account 
for the three-dimensional nature of the structure. 



plained recent experimental reports of features such 
as narrow-band resonances near the band edge' 3 ' 14 and 
showed excellent agreement with measurements on GaAs 
PC structures also presented. Similar theoretical findings 
were later reported and confirmed by Mazoyer et al. . 
In this work, we present and expand on the theory ex- 
ploited in Ref. 1 3 and provide a full derivation. Specifi- 
cally, we introduce a non-pcrturbative theory of coherent 
optical scattering over multiple periods of a disordered 
waveguide instance. The theory combines Green func- 
tion techniques and coupled mode formalisms with wave 
amplitudes calculated at each point along the length of 
the waveguide, where coupling coefficients include the 
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full three-dimensional disordered structure. In Section I, 
we introduce the theoretical formalism and derive the 
coupled mode equations for the forward and backward 
propagating Bloch fields. In Section II, we discuss the 
disorder model, and Section III implements the model 
with examples of simulated PC waveguide transmission 
and forward wave intensity. Finally, we conclude in Sec- 
tion IV. 



I. THEORY 

A. Waveguide Bloch Modes 

The ideal PC waveguide is periodic along the propa- 
gation direction (x) with periodicity a: e(r + ax) = e(r), 
where e(r) is the dielectric constant that we will assume 
is real and x is a unit vector. Consequently, Bloch's The- 
orem applies and the electric field mode may be written 
as Efe (r) cx e k (r) e lkx , where k is the Bloch wave vec- 
tor and efc(r) is the periodic Bloch mode. The magnetic 
Bloch mode (r) is defined similarly. Due to the Hcrmi- 
tian property of the Maxwell wave equations, the Bloch 
modes are orthogonal and, using the electric field modes, 
can be normalized through 16 
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where 5 k ,k' is the Kronecker delta; a similar relation holds 
for the magnetic field. The use of this relation as a pro- 
jection operator requires integration over the volume of a 
unit cell. For the present work, since we are interested in 
developing sub unit-cell propagation equations, we would 
prefer the integration was over only the plane perpendic- 
ular to the propagation direction. Using the electric and 
magnetic field orthogonality relations, the Maxwell con- 
stitutive relations, and the divergence theorem, one can 
derive 
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where the integration here is performed over a single 
plane transverse to the propagation direction. For k ^ k', 
the term in brackets is non-zero and the integral must 
evaluate to zero. For k = k' , the integral can be recog- 
nized as the power flux at the transverse plane which is 
clearly non-zero (except for a radiation mode propagat- 
ing perpendicular to the slab). Thus, a new projection 
(orthogonality) operator can be defined as 17 
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where E p (r) is the field being projected and x = xo is 
an arbitrary plane. This result is in agreement with 



that of Marcuse 18 and the standard form for overlap 
integrals 1 ' 1 . The projection operator Vk has the useful 
property that V k e k >(r)e lk ' x = 5 k ,k>- 



B. Green Function Approach for the Electric Field 

The electric-field properties of the disordered structure 
can be calculated analytically from Green function solu- 
tion to the electric field wave equation, namely 



E(r;w) =Ei(r;w) 



dr'G(r,r'; W ).^^, (4) 
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where P(r';w) is the polarization density due to the dis- 
order in the system (defined later), Ej(r; oj) is the electric 
field in the ideal system, and G(r,r';w) is the photon 
Green function where the overbar represents a tensor or 
dyadic. The Green function is a dipole solution to the 
Maxwell wave equation: 

V x V x - e (r) G(r,r';w) = 5(r-r')T, 

J5) 

where 1 is the unit dyadic. For convenience, we partition 
the Green function into contributions from the bound 
waveguide mode, radiation modes, and other modes as 

G(r, r'; w) = G B (r, r'; w) + G R (r, r'; w) + G G (r, r'; w). 

(6) 

The bound mode Green function is given analytically 
from properties of the bound mode' 1 ' 17 
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where the group velocity, v g , is assumed positive (in the 
case of anomalous dispersion, k is then negative), ® is a 
tensor product, e_fe(r) = e^(r), and 9(x) is the Hcav- 
iside step function, equal to 1 if x > and is x < 0. 
The mode properties can be calculated with any mode 
solving technique; for example, we use a freely available 
plane wave expansion code 20 . 

The radiation Green function, Gj{(r, r;u), contains 
contributions from the continuum of radiation modes 
above the light line that are not confined to the slab by 
total internal reflection. The radiation Green function, 
whose contribution is significantly smaller than the dom- 
inant bound mode, is rather featureless and is well ap- 
proximated by using a homogeneous dielectric slab with 
an effective permittivity determined through numerical 
FDTD simulations. We compute the radiation Green 
function efficiently by using the method of Paulus et al. 
(see also Ref. 17 for more details of our specific imple- 
mentation) . 

The remainder of the contributions to the Green func- 
tion are contained in Go(r, r'; uj) ('O' represents others), 



3 



such as the possibility of having other modes (bound or 
leaky), and the divergence contribution of the real part of 
the Green function as r — > r'. Since we consider a wave- 
guide with one bound mode in the frequency range of 
interest, we can safely neglect other bound modes. For 
the divergent contribution to Go(r,r';w), we shall ne- 
glect its contribution in this work; the dominant effect 
is to cause a ridged frequency shift"" and introduce local 
field corrections 23 ' 24 . 



C. Forward Wave Envelope Equation 

The electric field in the ideal waveguide can be decom- 
posed into the complete Bloch-mode basis consisting of 
the target bound waveguide modes e±fc(r), and the set 
of radiation modes {q(r)} as 



forward wave becomes 



E(r; u) = So [e fc (r) e ikx Vf (x) + ej(r) ■ 



where So is an amplitude and ijjf(x), V'b(x), and {ipq(x)} 
are the envelopes for the forward, backward, and ra- 
diation modes. We stress that we use envelopes only 
for convenience and do not require that they are slowly 
varying. We are only interested in the envelopes for the 
bound waveguide modes but we initially track the radia- 
tion modes to include radiation scattering. 

The field in a disordered waveguide can be calculated 
analytically from Equation 4, using the effective PC 
waveguide Green function and the disorder polarization 
density P(r;cj) = £ Ae(r) E(r; w), as 

E(r;w) ~E ?; (r;w) 

+ y"dr' [G B (r,r'; W ) + G R (r,r';w)] 
•[A £ (r')E(r';a;)], (9) 

where Ae(r) = e(r) — £i(r) is the disorder function and 
£i(r) is the dielectric constant for the ideal structure. We 
assume an initial electric field Ei(r;w) = So e k (r) e lkx , 
and a total field including scattering E(r;w) given by 
Equation 8. 

We begin by projecting Equation 9 onto a forward 
propagating wave by operating with Vk ■ We then multi- 
ply by Sq 1 and differentiate with respect to x. The left 
hand side becomes simply d^/>f(x)/dx. The projection of 
Ej(r; lu) equals 1 and differentiating eliminates the contri- 
bution of the field in the ideal structure. This derivation 
will transform the integral description of the total elec- 
tric field into a set of coupled propagation equations and 
the electric field in the ideal structure will be included as 
a wave injected from the input port. Equation 9 for the 
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The terms on the right hand side all arise from the 
projection of the Gs(r, t';u>) term; the projection of 
the Gn(r,r';w) term is since the constituent radia- 
tion modes are orthogonal to the chosen bound mode. 
The volume integral has been converted to an integral 
over the transverse plane by the derivative of the Heavi- 
side function in Gb(i", r';w). The scattering coefficients, 
corresponding to forward-forward, forward-backward, and 
forward- radiation scatter, are 
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An analogous equation to Equation 10 for dV'b(x) / dx is 
formed by projecting Equation 9 onto a backward prop- 
agating wave. One has 
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where the negative sign arises from the Heaviside func- 
tion in Equation 7, Cb b (x) = Cff(x), Cbf(x) = cj^(x), and 



Cbq(x) = — / / dydze k (r) e lkx ■ q(r) e lk « x Ae(r). 



D. Disorder-Mediated Coupled Mode Equations 

Next, we seek to eliminate the ip q (x) from the equa- 
tion since there are a large (infinite) number of radiation 
modes, and we would rather not have to solve for all the 
V'q(x). We project Equation 9 onto any one of the radi- 
ation modes to derive a radiation mode envelope equa- 
tion. The left hand side becomes simply ipqix). Only 
the Gn(r,r';w) term on the right hand side will have a 
non-zero projection since any chosen radiation mode will 
be orthogonal to the bound waveguide modes. Thus we 
obtain a set of equations, one for each of the radiation 
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modes q, 
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(15b) 
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(15c) 

There are three sources of energy for the radiation modes 
that are expressed as three terms on the right hand side 
of Equation 15: scattering from the forward wave (15a), 
scattering from the backward wave (15b), and scattering 
from all the radiation modes (including self-scattering 
from the current radiation mode into itself) (15c). 

First, we omit 15c, since we assume that scattering 
is just a loss mechanism and inter-radiation-mode scat- 
tering will not feed back into the waveguide modes. 
We also neglect 15b; this would give rise to radiation- 
assisted back-scattering where light from the backward 
mode scatters into a radiation mode and then the for- 
ward mode. These assumptions are reasonable because 
the radiation modes quickly leak from the slab and so 
do not interact with the scattering regions for very long. 
This leaves only 15a which accounts for loss from the for- 
ward mode into the radiation modes. The Vq prefix in 
Equation 15 is a projection operator acting on the radia- 
tion Green function. In Equation 10, the projected Green 
function (in ipq) is multiplied by the basis vector (in Cf q ). 
Since the set {q(r)} spans all radiation modes included in 
G R (r, r';w), this is an identity transform of G R (r, r';u>) 
and Equation 10, under substitution by Equation 15, be- 



v g —i/)t(x) =ic s (x) ipt(x) +icfb(x) e l2kx ipb(x) 

+ icf r (x)ipi(x), (16) 

where the radiation coupling coefficient Cf r is given in 
Equation 13 (which is further simplified below). Note 
that we have conveniently eliminated the sum over q. 

For the backward wave, Equation 14 is transformed 
using Equation 15 with only term 15b retained. The 
backward wave equation is 
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The final coupled mode equations are Equations 16 
and 17. The coupling coefficients can be physically in- 
terpreted as eg = Cbb (11) driving scattering from a 
mode into itself, Cbf = cj^ (12) driving scattering into 
the counter-propagating mode, and Cf r and Cbr driving 



scattering from the waveguide mode into radiation modes 
above the light line. With the elimination of the radiation 
mode envelopes, the coupling coefficients into radiation 
modes (e.g., 13) become 
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x e ikx e fc (r) • G R (r, r'; uj) ■ e* k (r') e - lkx . (19) 



Importantly, this theory incorporates the full three- 
dimensional structure of the waveguide, Bloch modes, 
and disorder functions in calculating the scattering. 

The radiation scattering coefficients of Equations 18- 
19 are difficult to evaluate due to the integral over the 
entire waveguide. Although we assume disorder between 
holes is uncorrclatcd in the expectation sense, for any in- 
stance of disorder, there may be a non-zero correlation 
between holes mediated by radiation modes. However, 
we are primarily interested in coherent scattering that is 
contained within the waveguide, and can reasonably as- 
sume that any field scattered out of a bound mode will 
not be scattered back into a bound mode; this is justified 
as the bound mode scattering channel is by far the domi- 
nant one. Therefore, we can simply the radiation loss by 
using Cf r = i (a ra d) v g /2a where (a r ad) is the incoherent 
average radiation loss 6 
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Comparing Equations 20 and 13, the former is just the 
expectation value of the imaginary part the later inte- 
grated over a unit cell. The factor of 2 is necessary to 
convert from a power loss to an amplitude loss. 

For modelling an incident field at one end of the wave- 
guide, the boundary conditions for a wave injected into 
the waveguide (and consistent with Ej(r;cj)) are 



^(^start) = 1, 
iM^end) = 0, 



(21) 
(22) 



where a; sta rt and a; en d are the positions of the input and 
output ports. The propagating envelopes are then com- 
puted at all spatial position within the waveguide using 
the presented coupled mode equations (Eqs. 16-17). Wc 
stress that the full three-dimensional Bloch mode and dis- 
ordered holes are self-consistently included in these final 
coupled-mode equations. 



II. DISORDER MODEL 

The equations can now be used with any disorder 
model. In our experience 4 ' 11 and in agreement with the 
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FIG. 2: a) Schematic of a hole with disordered perimeter 
and straight side walls. We describe the statistical properties 
of the disorder with a RMS roughness a and a correlation 
length l p measured around the circumference, b) Example of 
a disordered hole profile used in the calculation (blue). The 
ideal radius (dashed black) and correlation length (short red 
arc) are shown for reference, and R indicates the nominal 
radius of the unperturbed hole. 



an integration, where /(r^, fa) represents one of the fields 
and is slowly varying over the relevant length scale. The 
field f(ri,<fii) can be expanded in a Taylor series along 
the radial coordinate to evaluate the integral as 

Jdn Ae^fj, fa) f(n, fa) 

= Jdn Ae(n,fa) 

x (f(R, fa) + f(R, fa)( n -R) + 0((n ~ R) 2 )) 
= .f(R,fa) J d n Ae(n,fa) 

+ f'(R, fa) J d n As{n,fa)(n -R) + 0((r< - R) 2 ) 
= f(R,fa)(e 2 -e 1 )Ar i (fa) 

+ f'(R, fa) (£2 - ei) Ar ^ )2 + 0(An(fa) 3 ). 

(25) 

To include the disorder to first order in Ari(fa), it 
is sufficient to take the field at the ideal hole radius 
f(R,fa). For convenience of notation, we then rewrite 
Equation (24) as 

Ae t {n, fa) = (e a - £i) S(n - R) A n {fa), (26) 

so that 

Jdn Asi(ri,fa) f{n, fa) = f{R,fa) O2 - £1) Ar^fa), 
which agrees with Equation (25) to first order. 



analysis of images of PC slabs 25 , we have found that dis- 
order in PC slab structures is dominated by perturba- 
tions of the perimeter of the holes, as shown in Figure 2. 
We take the radial perturbation Ar to be a Gaussian ran- 
dom variable with a mean of and a standard deviation 
of a. Two radial perturbations are correlated by 

(An{fa)Ar 3 {fa } )) = a 2 e -*l**-*}l/'p 6 hJ , (23) 

where the subscript indexes the holes, fa is the angular 
position of the point measured about the centre of the 
hole, R is the ideal hole radius, and l p is the correlation 
length measured around the circumference. 

The change in dielectric constant about a single hole i 
is given exactly by 

Ae i (r i ,fa) = (e 2 -e 1 ) [0(r< - R) 

- G(n -R- An (fa))], (24) 

where (ri,fa) are cylindrical coordinates centred about 
hole i. This form holds for both positive and negative val- 
ues of Arj(fa). The disorder Ae appears in the formalism 
in spatial integrals where it is multiplied by functions of 
the electric fields and Green function. We consider such 



III. IMPLEMENTATION 
A. Ideal Structure 

This calculation requires, as inputs, the ideal wave- 
guide mode dispersion and spatial field distribution. As 
a representative example we consider a Wl semiconduc- 
tor waveguide with pitch a = 480 nm, slab thickness 
h = 160 nm, hole radius R = 95 nm, and index of re- 
fraction n = 3.18. The dispersion of the waveguide mode 
is shown in Figure 3a) (blue, solid, left scale) along with 
the group index (green, dashed, right scale). Near the 
band edge (k = 2ir/a), the group index is large, increas- 
ing scattering as the light slows down. The spatial dis- 
tribution of the electric field in the centre of the slab is 
shown in Figure 3b). 

B. Numerical Implementation 

To solve Equations 16 and 17 numerically, the cou- 
pling coefficients are assumed to be constant over a short 
(Ax <C a) interval in x and are integrated analytically. 



6 



a) 




y 



FIG. 3: Properties of the nominal structure, a) Dispersion of 
the waveguide mode (blue, solid, left scale). The continuum 
of radiation mode above the light line is indicated by the 
shading on the left side of the figure. The group index (green, 
dashed, right scale) is also shown and diverges at the band 
edge (k = 2n/a). b) Distribution of the transverse component 
of the electric field Bloch mode at the middle of the slab near 
the band edge. 



This yields a pair of transfer equations linking the en- 
velopes on either side of the chosen interval. In this way, 
a set of transfer equations that span the entire waveguide 
length can be built, and then solved using linear algebra 
techniques. This approach is particularly amenable to 
adding reflective facets and other features by simply in- 
cluding an appropriate transfer matrix. 

The average coupling constants for each interval are 
calculated by, for each hole, generating an instance of 
a disordered profile from the statistical distribution of 
Equation 23. The coupling coefficients are calculated at 
multiple points within the interval, and then averaged. 
Typically, there arc 20 intervals per unit cell to satisfy the 
assumption that the coefficients are relatively constant. 
As shown in Figure 4, if the discretization of the unit cell 
is too coarse, the loss is underestimated. Thus, one must 
include sub unit-cell propagation effects. 

We highlight that the calculation is orders of magni- 
tude more efficient than standard brute-force numerical 
techniques, e.g., FDTD. We also note that we only need 
to calculate the coupled mode coefficients wherever disor- 
der has an influence, namely at the hole interfaces. How- 
ever, the final computation, though efficient, is not in- 
stantaneous. Producing a high resolution transmission 
spectrum (1000 frequency points) for a 1 mm waveguide 
(2 500 unit cells and 50 000 grid points) takes approxi- 
mately 1 cpu day (on a 2.4 GHz AMD Opteron proces- 
sor). However the calculations at each frequency are in- 
dependent and the total calculation can also be greatly 
accelerated by exploiting parallelism. In contrast, we es- 
timate that a minimum of about 40 GBytes of memory 
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FIG. 4: Mean transmission through 500 disordered waveg- 
uides (blue, solid) as a function of the number of intervals 
each unit cell is divided into. The error in the mean is marked 
by the dash-dotted limits and the mean agrees well with the 
prediction of the incoherent calculation (red, dashed), except 
for very coarse discretizations. 

and 5800 cpu days are required to perform the simula- 
tion using FDTD. Clearly, this semi-analytic treatment 
is a significant advantage. 



IV. COMPUTED TRANSMISSION SPECTRA 

Figure 5 shows transmission spectra for four disordered 
waveguides calculated by solving Equations 16 and 17 
(blue, solid). For reference, previous incoherent scatter- 
ing results, computed within a second-order Born approx- 
imation ,, J are also shown (red, dashed); we also note that 
extensions to the incoherent scattering theory to account 
for multiple scattering have been introduced recently 11 . 
Each row of plots is for a different waveguide with the left 
plot showing a broad frequency range and the right plot 
showing a narrow frequency range near the band edge. 
The top row is for a disordered 1.5 mm waveguide (3125 
unit cells). The second row is for a different disorder 
instance of the same 1.5 mm structure. Experimentally, 
this would be similar to carrying out measurements on 
a second waveguides fabricated with nominally identical 
parameters. It has the same qualitative shape but the 
particular disordered resonances are substantially differ- 
ent. This is important if it was desired to take advantage 
of these sharp resonances since their resonant frequency 
cannot be easily designed. The third and forth rows are 
for the same disorder instance as the second but with the 
length reduced to 1.0 mm and 0.5 mm respectively. Here 
the qualitative roll off changes due to the length reduc- 
tion but disordered resonances can be found at similar 
frequencies across the three lengths, especially between 
the 1.5 mm and 1.0 mm cases. 

We can examine the position-dependent distribution 
of energy in the waveguide under c.w. illumination. In 
the second row, right column of Figure 5, a neighbouring 
transmission minimum and maximum are marked with 
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FIG. 5: Simulated transmission spectra of four disordered Wl 
waveguides using the new coherent scattering theory (blue, 
solid) and the first- and second-order Born incoherent theory 
(red, dashed)". Each row of plots is for a different waveguide 
with the left plot showing a broad frequency range and the 
right plot showing a narrow frequency range near the band 
edge. Plots (a) and (b) are for a disordered 1.5 mm waveguide. 
Plots (c) and (d) are for a different disorder instance of the 
same 1.5 mm waveguide. Plots (e-f ) and Plots (g-h) are for the 
same disorder instance as (c-d) but with the length reduced to 
1.0 mm and 0.5 mm respectively. The calculation uses a RMS 
roughness of a — 3nm, and a disorder correlation length of 
l p = 40 nm. The forward wave intensity as a function of 
position is given in Figure 6 for the two points marked with 
crosses in d). 



red crosses. The forward wave intensity at these frequen- 
cies is plotted in Figure 6. Although the points are very 
close in frequency, the minute difference in group index 
(ri g = 25.11 compared to n g = 24.96) creates a difference 
in the accumulated phase and a dramatic change in the 
transmission. 

By including multiple, coherent scattering we repro- 
duce the experimental phenomenon of sharp spectral 
resonances near the band edge. Although initially 
unexpected, these features are just Fabry-Perot-likc 
fringes between extrinsic scattering sites. The slow 
group velocity enhances scattering to create the 
scattering sites and also increases the effective cavity 
length between sites, narrowing the resonance line- width. 
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FIG. 6: Forward wave intensity in a disordered waveguide at 
two wave vectors. The blue curve (n B = 24.96) corresponds 
to a local transmission maximum and the green curve (n g = 
25.11) is a neighbouring transmission minimum. These two 
curves correspond to the red crosses in the Figure 5d). 



V. CONCLUSIONS 

We have described and applied a theory for self- 
consistcntly modelling coherent scattering in a disordered 
PC waveguide instance, allowing one to map directly onto 
a realistic experimental situation. Slow light propagation 
enhances back scattering (and, to a lesser extend, radi- 
ation scattering) leading to high losses near the band 
edge. The formation of sharp spectral resonances near 
the band edge is shown which is mediated by Fabry- 
Perot-like resonances between disorder sites. This the- 
ory is computationally efficient, making the analysis of 
very long waveguides (thousands of periods using the full 
three-dimensional structure) feasible on a desktop com- 
puter. Although the presented model may not be quan- 
titatively exact (e.g., it neglects local field effects), the 
qualitative results such as the formation of sharp reso- 
nances near the band edge certainly can, and already 
have been, used to explain a rich range of experimen- 
tal features without introducing any fitting parameters 1 ! . 
The role of local field effects will be reported in future 
work, and the effects on incoherent frequency shifts are 
described elsewhere 26 . 
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